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The behavior of two interacting populations, "hosts" and "parasites", is investigated on Cayley 
trees and scale-free networks. In the former case analytical and numerical arguments elucidate 
a phase diagram for the Susceptible-Infected-Susceptible model, whose most interesting feature is 
the absence of a tri-critical point as a function of the two independent spreading parameters. For 
scale-free graphs, the parasite population can be described effectively by its dynamics in a host back- 
ground. This is shown both by considering the appropriate dynamical equations and by numerical 
simulations on Barabasi-Albert networks with the major implication that in the thermodynamic 
limit the critical parasite spreading parameter vanishes. Some implications and generalizations are 
discussed. 
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I. INTRODUCTION 

Population models, or reaction-difFusion systems have 
attracted enormous interest both in the statistical physics 
community and as abstract versions of real biological dy- 
namics. One particular aspect is the presence of phase 
transitions and the contact process or directed percola- 
tion in various disguises (see below, 0,0). 

Host-parasite or predator-prey systems are a natural 
extension of single species models. By their classical re- 
sults Lotka and Volterra were able to explain the na- 
ture of abundance oscillations in interacting species [^0. 
In regular landscapes or lattices, with a finite spreading 
rate of the species, these oscillations appear as travel- 
ing waves, which can be regular or chaotic, depending on 
the interplay of time scales in population dynamics and 
spreading, though it is not clear if thephenomenon sur- 
vives m the thermodynamic hmit [1 S S H H E Ell ■ 
In nature they have been observed in different systems, 
to name two extreme cases, e.g. in vole populations ^3 
and for human diseases such as measles [l^ . In the case 
of measles in a population living on a landscape of non- 
trivial island structure, power law fluctuations are found 
instead [l^. 

Much of these ideas have recently been generalized in 
the context of small-world or in particular "scale-free" 
graphs El El El d. For the latter, a perfectly valid 
example is given by epidemics of viruses in the Internet 
since it has as a graph a fat-tailed probability distribu- 
tion of the number of nearest neighbors, P{k). Recently, 
various models have been studied as the particulars of the 
structure — as the so-called degree distribution gamma 
in P{k) — are varied. A fundamental discovery 

concerning disease spreading is an absence of epidemic 
threshold in the limit of infinite graphs and the finite- 
size effective "critical point" obero an unusual scaling as 
L, the graph size, is varied |19ll20ll21j . 

This closely relates to the present work where we study 



the influence of a network or graph like structure of the 
underlying landscape on host-parasite or predator-prey 
dynamics. The main findings are (i) the absence of os- 
cillations, (ii) the absence of an infection threshold in 
the limit of an infinite scale free graph, and (iii) the ex- 
istence of two separate transitions in the case of Bethe 
lattices with finite coordination number z ("empty" — > 
"hosts only" , "hosts only" — > "hosts plus parasites" , but 
no transition "empty" — > "hosts plus parasites"). The 
structure of the rest of the paper is as follows: Section II 
contains the necessary definitions, and the two following 
ones analytical considerations and, to compare, numeri- 
cal simulations of the models. Finally, Section V finishes 
the paper with a discussion. 



II. MODEL FORMULATION 

A. States and rates 

A basic model for epidemiological applications is the 
contact process, or the so-called SIS model. Here one 
considers individuals living on the nodes of an underly- 
ing graph which are either infected (/) or susceptible {S) 
to an infection. An infected individual may spread the 
disease to a susceptible one if both are in contact, i.e. 
if they live on neighboring nodes of the graph. Infected 
individuals recover with a certain rate and in this simple 
version immediately become susceptible for a new infec- 
tion. So the dynamics of the SIS model is defined by the 
rates 
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if any neighbor is infected 



(1) 



In this work we generalize the SIS model to a system 
with hosts and parasites (HP). In other words we consider 
infections of a second kind only able to spread onto sites 
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FIG. 1: States and rates. 

with infections of first kind. So each node in the graph 
can be in three possible states: Empty (e), populated 
by a healthy host (h), or by a host with parasite (p). 
Between three states there are six possible transitions so 
the dynamics are defined by the following rates 

r^^fi = A if any neighbor has a host (h) 

rfi—tp — a II if any neighbor is parasitized (p) 

rh-.e = 1 (2) 

re^p = 

Tp^h = e 

— >e — 

As in the SIS model defined above the decay of the host 
or first kind infection sets the time scale {vh^e = !)■ In 
biological systems a > 1 (even ^ 1) if the parasite affects 
the health of the host. A benefit would mean a < 1. We 
shall consider cases in which the parasite virtually does 
not die "on its own" but only when the host is killed, i.e. 
the case « e ^ 1, A, a/i, a. 

In Section [nil we present approximate analytical solu- 
tions following [13 to the the model of Eqs. (|3J) which 
are compared to Monte Carlo simulations in Section ITVI 
Particular interest lies in parasite extinction and its de- 
pendence on the parasite spreading rate a/x. But first we 
define the types of graphs used in our simulations and 
calculations. 



B. Graphs 

We study the population dynamics of the HP-model 
on two types of graphs, on Bethe lattices and on scale 
free Barabasi- Albert (BA) graphs, in their standard ver- 
sion [231 . A Bethe lattice of coordination number z is 
an infinite tree, where each node has z neighbors. When 
constructing a finite lattice, or Cayley tree, starting from 
a central node with z neighbors and adding 2 — 1 new 
neighbors to each boundary node, the number of bound- 
ary nodes grows exponentially. It therefore remains a 
finite fraction of the total number of nodes in the finite 
tree, which makes this construction unsuitable for Monte 
Carlo simulations. 

This difficulty can be overcome by a slight modifica- 
tion 22] where a sparse homogeneous graph that closely 
approximates the Bethe lattice without any boundary 



nodes is constructed. Take L nodes and label them by 
integers from 1 to L. Connect node i to node [i + l) for 
each i and connect node I to node L. Construct {z — 2) 
independent random pairings of the nodes (an easy way 
to construct pairings is to sort the nodes randomly and 
pair the first node of this new order with the second one 
etc.) of the nodes and add an edge for each pair. By this 
procedure, we get a graph in which each node is of de gree 
z. For large enough graphs, the loops are negligible |23| 
and this is a sufficient approximation of a Bethe lattice. 

Here, we also use the standard version of Barabasi- 
Albert graphs (BA-graphs) [23|. These are constructed 
as usual. New nodes are added one by one connecting 
them with m = 3 links to the previous ones. From these, 
the neighbors are chosen with a probability proportional 
to their respective number of links (preferential attach- 
ment) . By this construction highly linked nodes are likely 
to obtain even more neighbors as the graph grows, which 
results in a "fat tail" distribution of probabilities for a 
node to have coordination number fc, P{k) ^ pH . 
The BA-graphs have very weak degree correlations, i.e. 
the conditional probability for a node of degree k to have 
a neighbor with k' is rather trivial [TEj compared to many 
other models and real networks. 



III. MEAN FIELD AND DOUBLET 
APPROXIMATION 

A. Bethe lattice 

1. Singlet (mean field) approach 

In this subsection we extend the known solution for 
the SIS model on a Bethe lattice to the HP model. 
Ph and pp denote the density of hosts and parasites, re- 
spectively. For simplicity we consider the limit e = 0, 
so parasitized patches do not supply host individuals to 
neighboring empty patches. The rate equations for the 
densities can be written as 

dtPh = -Ph + >^ii-Pp-Ph) Q~apph^ 

dtPp = -a pp + ap Ph <& (3) 

with $ = 1 - (1-pp)^ and e = 1 - (l-p^^)^ 

In the absence of parasites the host population follows 
the dynamics of a SIS model. The trivial state ph = 
is stable for A < 1/z and unstable otherwise. In other 
words, the host population can survive only for A > 1/z. 

Similarly, the pure host phase is stable if parasites can- 
not invade, i.e., if the growth rate of a small parasite 
population is smaller than its death rate, 

a p z Ph < a or p < ^— = u'^"'. (4) 

zph 

From this formula it can be seen that p"^'^ — > 00 as ph 
0, so there is no "tricritical point" in the phase plane, 
beyond which a direct transition from the absorbed state 
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to the coexistence state can be seen. The phase diagram 
is drawn in Fig. |21 



2. Doublet approach 

The singlet approach neglects occupancy correlations 
between adjacent nodes. The next logical step is a pair 
or doublet approximation which explicitly treats the joint 
probabilities to find two unparasitized hosts next to each 
other {Phh)i a healthy host next to a parasitized one 
{Php), and two parasitized next to each other (Ppp) in 
addition to ph and pp. This approximation is widely 
used, we want to emphasize its application to a spatially 
uniform insect host-parasitoid model 23, '24], to the con- 
tact process in a one-dimensional chain ^ and in general 
over a wide class of models . 

The approximation uses the probabilities Paa' to find 
the nodes adjacent to a randomly picked bond in states a 
and a' G {e, h,p}, as well as the conditional probabilities 
p^rjjr' to find a randomly chosen nearest neighbor of a 
cr'-node in state a. Three-point and higher correlations 
are neglected, so the conditional probabilities to find a 
(7-node next to a cr'-node which is itself linked to a third 
node with state a" are approximated by 



Pa\a'cr" ~ Pa\a' V Cr". 

From there one obtains the rate equations 



dtPh 
dtPp 
dtPhh 



-1 - zappp\h + z\p^\h 



Ph 



-a + zapp^p 



Pp 



2 + 2(z-l) ap pp\h 



Phh 



+ A 



1 + {z~l)ph\. 



Ph 



(5) 

(6) 
(7) 

(8) 



dtPpp 



dtPhe 



dtPpe 



dtPee = 



dtPhp = -jl + a + a/i 1 + (z-l) ppih ^ Php (9) 

+ (z-1) \ph\ePpe + 2(^;-l) apPp\hPhh 

-2aPpp + ap 1 + {z-l)pp\h Pi 
-|l + {z-l)appp\h 



hp 



(10) 



+A 



+ 2(Z-1) A Ph\ePee + "^Phh + aPhp 



P 



pe 



- a + {z-l)Xph\^ 
+ {z-l)ap Pp\hPhe + Php + 2aPpp 

-2{z-l)Xph\ePee + Pie + aPpe 



(11) 

(12) 
(13) 



The joint probabilities Pa-a-' can be expressed in terms of 
the conditional probabilities as 



Paa' = PaPa'\a{^ - 



(14) 



where ^o-,ct' is the Kronecker symbol. The factor 2 for 
a ^ a' refiects the two possible choices, because a can 
be on either end of the bond. 



There are some subtleties in Eqs. (jH)l- (|13|) that might 
not be immediately obvious. In Eq. ||SJ), for instance, 
there is a factor of 2 in the first term. That term describes 
a process where an edge connecting two host-sites turns 
due to a death of a host into an edge connecting an empty 
site to a host-carrying site. The prefactor comes from the 
fact that this can happen in two ways, i.e. either of the 
the two hosts can die. For similar reasons, a prefactor 
of 2 can also be found in the second term of Eq. ((SJ. 
However, the rest of the terms in the equation do not 
have these prefactors since a similar symmetry does not 
exist. 

In principle Eqs. l|B|l- H13() are solvable in the steady 
state. Consider first the SIS model, i.e. the case without 
any parasites. Setting p = Q and looking at the steady- 
state of Eq. lO immediately yields 



Pe\h 



(15) 



Similarly, setting Ppe = in Eq. (|13|l and using the iden- 
tities Pee = PePe\e and Phe = 2pePh\e givCS 



Pe\e 



1 



(z-l)A 



Expressing ph as 



Ph = Pe = (1 - Ph) 

Pe\h Pe\h 



(16) 



(17) 



using the identity Pe\e + Ph\e = 1; and plugging in 
Eqs. (|15|l and (|15|) finally gives 



Ph 



(z-l)A-l 
(z-l)A-l/z' 



(18) 



from the numerator of which the critical point follows: 

1 



A° = 



(z-1) 



(19) 



Note that this is different from the mean field result 
^MF = It is also worth noting that rigorous math- 
ematical results of the contact process j2^ give bounds 
on the critical point as 



1 ^ 1 

- < Ae < 

z z — 



1 



(20) 



Next consider the boundary between the parasite- 
absorbing and the coexistence phases. Here, hosts live 
well while parasites are near extinction. Expanding the 
steady state solution in the limit of small parasite popu- 
lation we derive an equation for the phase boundary: De- 
fine two auxiliary quantities A — p^p and B — Ph\p+Pp\p- 
Form an equation for dtA and set it to vanish since we 
are looking at the steady-state 



dA 
'dt 



d_Php 
dt 2pp 



1 

2pp 



dP, 



hp 



dt 



PhpMp = (21) 
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where the rate equation (O has been used, and Mp is the 
Malthusian parameter or growth rate at low densities of 
the parasites, i.e. 



afizA . 



(22) 



Plugging Eq. ® into Eq. using Eq. (O and the 

results of Eqs. lfT3|l and (fTB|) at vanishing parasite popu- 
lation we arrive at 



2z{z-l){l- B)X^ + 
+{2z^afiA{l - A)+ 
-2z[B -I- A{1 + 2a^i)]}\- 

-2(z- l)Aa^i = 0. 



Similarly, starting at the rate equation for B, 

9pe\p 



dB 
'dt 



■VP 1 / dPe„ 

}£. — I P M — 

dt " 2p„ V'"^^^" dt 



(23) 



(24) 



using Eqs. 1)12(1 and Eq. H14|) together with the results of 
Eqs. ((T5)l and ifTBI) . one gets 

2z{z-l){l- B)X^ + 
+2z'^anA(l - B)X+ 
+2z[{B - A){1 - a) - 1]X- 

-2{z-l)Aap = (25) 

given that X ^ 0. 

Now, solve for A in the steady-state version of Eq. ||7J), 
substitute this to Eqs. (|23|l and (|25(l and eliminate B from 
the resulting two equations to get 



z{z — 1)X^ + zaX 



z{z - IfX^ + z[{z - \){a - 1) - a\X + a(l - z) 

(26) 

for the phase boundary between parasite-absorbing and 
co-existence phases in the (A, /^)-plane. Note that, con- 
trary to the mean-field approximation, the phase bound- 
ary defined by this equation does not meet that defined 
by Eq. H19I) at /i ^ oo since A — X^ is not a zero of the 
denominator of Eq. H2()|) . It also holds that p,c — >■ ^/{z — l) 
as A — > cxD so that in this limit the parasites always find 
hosts on all nodes and therefore behave as the SIS model 
does. 

In addition to the solution above we linearize the 
doublet rate equations around the previously obtained 
fixed point with a host population and no parasites, i.e. 



tional probabilities Po-'io- by joint probabilities 
in Eq. (|14|l we get a matrix which is of the form 



Ph 



p.. 



0. 



Replacing the condi- 
as 



M 



Mh M, 



hp 



Mr, 



(27) 



where Mh governs the stability of the "host only" solu- 
tion, Mfip the effect of a small parasite population on the 
hosts, and Mp the growth of parasites at low densities. 



The block in the lower left corner is zero since the state 
without parasites is an absorbing one, i.e. a perturbation 
in the host density cannot reintroduce parasite popula- 
tion into the system. 

The eigenvalues of a matrix with the structure of 
Eq. (|27|l are just those of Mh and Mp, irrespective of 
Mhp- The stability of the host population has been dis- 
cussed above, so we are only interested in the (real parts 
of the) eigenvalues of the matrix in the following equa- 
tion. 




-a 



2 

B 




A 

ap -2a 
V C 2a -a-X ) 



Pp 
Php 
P, 



.P 



pp 

ep i 



with 



B = (z-l)a^(zA-l)/(zA) - 1- 
C = (z — l)Q:/i/(zA) + 1 and 
A ^ (z-l)A-l. 



-a — a/i 



(28) 

(29) 
(30) 
(31) 



Note that A is proportional to the excess over the crit- 
ical host growth rate, A — Ac. 

The left column of the matrix in Eq. H28|) is empty 
except for the diagonal element, which gives the first 
eigenvalue, —a. We therefore restrict ourselves to the 
remaining 3x3 matrix. It is straightforward to calculate 
its eigenvalues explicitly, from which the phase boundary 
can be deduced as follows. For each fixed A, wc consider 
the real part of the largest eigenvalue of the 3x3 matrix 
as a function of and find its zero numerically, lead- 
ing to a point /i(A) that lies at the phase boundary. The 
results are shown in Fig.[21for the case a = 1.2 and z = 4. 

The absence of the tricritical point can be seen easily: 
As A \ Ac the excess growth rate A \ 0, and the matrix 
becomes lower triangular. All three diagonal elements 
yield negative eigenvalues, in particular in this limit B — > 
—apuj z — a — 1 < 0. In particular none of the eigenvalues 
approaches zero as — s- oo, which leads again to the 
conclusion that the two phase boundaries do not meet at 
this limit. 

In comparison to these results the mean field approxi- 
mation underestimates the critical values for the spread- 
ing parameters. It does not take into account the cluster- 
ing of populations, i.e., the fact that next to a populated 
site there is likely another one, which can not be invaded 
any more. So the possibility for growth is overestimated. 

The phase diagram of the HP model in the (^, A)- 
plane obtained from both theoretical approaches and 
from a stochastic simulation using graph approximation 
discussed in Sec. Ill Bl is drawn in Fig. [3 In the simu- 
lations, rough estimates for the phase boundaries were 
obtained by performing a series of simulations with dif- 
ferent A for each fixed [i and observing when the pop- 
ulation died out. The largest value of A at which the 
population dies out is then defined to be the estimate for 
the position of the phase boundary. From the figure we 
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FIG. 2: Phase diagram for Bethe lattice with z = 4 in 
the (A,/i)-plane with a = 1.2 and = 40000 nodes. Sin- 
glet approach (mean field) compared to doublet approach 
(pair approximation) both analytically and via stability anal- 
ysis and Monte Carlo simulations. The abbreviations denote 
the parasite-absorbing (p. a.), co-existence (c.e.) and empty 
phases, the last of which is the phase where both popula- 
tions will eventually become extinct. All three solutions are 
in qualitative agreement with each other. As one expects, the 
pair approximation predicts the need of higher growth rates 
(A and fi) than the singlet approach. 

see that both analytical solutions are in qualitative agree- 
ment with each other and with the numerical results. A 
property worth noting of the phase diagram is the lack 
of a "tricritical point" and thus the phase boundary be- 
tween empty and co-existence phases. Consider also that 
the singlet approach does reproduce the features of the 
phase diagram in the Bethe lattice case. 

B. Scalefree graph 

1. Singlet approach 

On graphs with non-constant degrees the occupancy of 
a node depends on its coordination number. In general, 
the higher the degree of a node, the greater is its tendency 
to be populated. Following Ref. [I^ the rate equations 
for the occupancies pjj and on nodes of degree k can 
be written as 

-pakptmk{X,p) (32) 
dtp'pit) = ~ap';{t)+fiakpUmk{\,p) (33) 

where 0/c(A,/i) and <l>fe(A,/i) are the probabilities that 
a given link points to an infected or a parasitized node, 
respectively. In Eq. (|32l) the first term on the RHS cor- 
responds to the death of the hosts, the second one to the 
host spreading and the third one to parasite spreading, 
diminishing the number of sites that carry host but no 



parasite. In Eq. (|33|l the first term on the RHS describes 
the death of the parasites while the second one encom- 
passes the spreading. It is known that there is no 
epidemic threshold if the distribution of node degrees is 
fat-tailed. 

The critical behavior of the HP model as obtained from 
the mean field equations above turns out to be incorrect 
and is in contradiction to the numerical findings. To see 
this, consider the rate equations (|32|l and l|33() in the limit 
of small p, i.e. by a Taylor expansion in p. The interaction 
term p,ap'^^k{X, P-) is quadratic in p since $fc(A, /i) ^ p 
and drops out from the expansion to first order. This, in 
turn, means that in this limit the host population behaves 
as in the SIS model and the parasite population dies out 
since its equation only has exponentially decaying solu- 
tions. Furthermore, this rules out the possibility of a 
zero epidemic threshold for the parasites, since when the 
spreading rate approaches zero also the prevalence does 
so. This leads to the aforementioned contradiction. The 
corresponding numerical results are presented in Fig. ^ 

2. Singlet approach with a substantial host population 

Next, we use the singlet approach to look at the be- 
havior of the parasites when the host population is well 
established. The calculation presented here is a straight- 
forward generalization of that in Ref. 0| . 

The rate equation of the parasites in a Markovian cor- 
related graph in the singlet approach can be written in 
the limit of small prevalence as 

d 

-glPp = -Pp + appUk , (34) 

where (f)k = J2k' kAkk'Pp , and Akk' = P{k'\k), i.e. the 
conditional probability that starting from a node of de- 
gree k and following a random edge one is lead to a 
node of degree k' . For uncorrelated networks, P{k'\k) — 
kP{k')/ (k) , where P{k) is the degree distribution of the 
underlying network. 

If the parameters are chosen such that there are plenty 
of hosts and that parasites are near extinction the feed- 
back coupling of the host population to the parasites can 
be neglected and pf^ approximated by a constant vector 
given by the solution of the SIS-model. The zero solu- 
tion Pp — Oyk is always a (formal) solution of the system, 
so we have to study its stability. Take Eq. (|34|l . denote 
P — (Pp ' ' ' Pp")'^ ^^'^ write the equation in a matrix form 

d 

— p = (-1 + afipf-^kAkk') P = {-I + apCkk')p , (35) 

where Ckk' = P^kAkk'- 

Looking at the matrix elements Ckk' gives 

P{k)_ P{k') 

^kk' T^ — i^k'k p— [ob) 

Ph Ph 
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where the detailed balance condition of the network l28 



kP{k'\k)P{k) = k'P{k\k')P{k') 



(37) 



has been used. From Eq. it follows that C and 
have the same eigenvalues since, if Vk is any eigenvector 
of C corresponding to eigenvalue A, then VkP{k)/ p{k) is 
an eigenvector of with the same eigenvalue. This, 
in turn, has the consequence that the spectrum of C is 
real. Again, the zero solution is unstable, if the matrix 
— / + ^.C has at least one positive eigenvalue, and the 
critical value of ^ is ^critical = A^"^ . 

Next, use the following corollary [23| of the Frobenius 
theorem. Let A]^y be any positive irreducible matrix. Its 
largest eigenvalue Am can be estimated from below as 



. / 1 

Tin y . 



Am > min { ^kk'i'ik') 



(38) 



where V'(^) is an arbitrary positive vector. Now, set 
i/;(fc) — kpf^ and A — C'^io get 



Efc- Y.iPlkPmpil P{k'\l)k'pt 

kpl 



Aa/ > min ■ 



mm ' 

k 



Y,ip\pmY.^'pi:p{k'\i) 



(39) 



Above /c„„ (Z, kc) denotes the average nearest neighbor 
degree of such neighbors that carry a host, conditioned 
that we are looking at a node of degree I. Since the aver- 
age nearest neighbor degree of all neighbors knn{l, kc) — 
^f,, k'P{k'\l) diverges [23 as fee — > oo and necessar- 
ily saturates to a constant value < 1 with large k, 

k\^l^{l, kc) must also diverge at the same limit. Thus the 
RHS of Eq. diverges, giving Aj\/ ^ oo and 



^critical 



at the thermodynamic limit. 



(40) 



connectivity k ik') is a {cr'), possible states being e, h or 
p. Q'^^i is the conditional probability that a randomly 
chosen edge that connects nodes with connectivities k 
and k' is such that the state of the node with connectiv- 
ity k' is <t' conditioned that the state of the node with 
connectivity k is cr. Let A^fc/ be as above. 

Using the notation above, the rate equations for the 
SIS model needed for the present treatment can be writ- 
ten as follows 



dtpt 



dtP^H 



-Ph 



\Y,k^kk'Pit 



r) pkk' 



xp; 



kk' 



(41) 
(42) 



-XY,^k'k"Pt7ik'-l)Q':', 



k 
h J 



where in Eq. (|42l) the first term on the right hand side 
denotes the process where an infected node gets cured, 
the second the process where a node of degree k infects a 
node of degree k' and the third the process where a node 
of degree k" infects a node of degree k', which in turn 
has another neighbor of degree k that is infected, turning 
the edge between the latter two into an edge connecting 
two infected nodes. 

For the HP model, only one rate equation is needed 
for the present treatment, namely that of the parasite 
prevalence 



dtPp = -app + fiaY^kAkk'Ptp 



(43) 



Now consider the steady-state in the SIS model. Mul- 
tiply Eq. H41|) by P{k) and sum over all k to get 



Ph = XPeh 



(44) 



Peh is the fraction of all edges in the network that connect 
an empty node to one with host and ph is the average host 
prevalence in the whole network. 

The last term on the right hand side of Eq. (|42|l is 
positive. Thus in the steady state we can write, leaving 
out the said term, 



P, 



kk' 



hh 



kk' 
eh 



(45) 



3. Doublet approach 



Multiplying this by kP(k)Akk' , summing over all k and 
k' and combining with Eq. we get 



Next, we formulate rate equations for a graph with a 
given degree distribution and degree-degree-correlations 
using the doublet approach or pair approximation. The 
correlations are included in the treatment since their use 
is natural in the context of pair approximations. The 
correlated network contains its uncorrelated counterpart 
as a special case. 

The notation is as follows: P^^, is the probability that 
a randomly chosen edge that connects nodes with con- 
nectivities k and k' is such that the state of the node with 



Phh > -^Ph , 

which implies for the relative density of host-host 
nearest-neighbor pairs that 



P, 



hh 



(Ph) 



1 1 

> - • > OO as Ph 

2 Ph 



0. 



(46) 



That is, in the limit of small population, the relative den- 
sity of host-host pairs is enormous. Thus the prevalence 
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correlations in nearest-neighbor nodes are also huge. 
Since the singlet approach neglects these correlations, 
this gives reasons to expect that it is not able to capture 
the properties in the HP model correctly, even though it 
is known that in SIS model it does [2?}. 

Consider Eq. H43(l in the steady state. Multiplying by 
P{k) and summing over all k gives 

Pp = pPhp (47) 

which in turn gives 

= — * oo as Pp , (48) 

PhPp f^Ph 

since /i — > as pp ^ 0. 

Eq. H47II tells that the number of edges through which 
the parasite population can spread is proportional to the 
parasite prevalence (instead of the product of parasite 
and host prevalences). This, in turn, tells that the dy- 
namics of the parasites is similar to the dynamics of the 
hosts in the SIS model (since in the SIS model the number 
of edges that can spread the population is proportional 
to the population density in the steady-state) and serves 
as an explanation to the zero threshold of the parasites. 

IV. MONTE CARLO SIMULATIONS 

For a numerical comparison we have simulated the 
host-parasite-model in Barabasi- Albert networks of sizes 
L = 2^^ . . . 2^^ under the conditions in which ph « 0.30 
and Pp <S Ph, i.e. with a stable hosts population and 
parasites close to extinction. The simulations are always 
started with random initial conditions by giving 25% of 
the nodes the status host and 5% of the nodes the sta- 
tus parasitized independently. Then the simulation is 
run for a given saturation period of 1000 MC-steps dur- 
ing which even the largest system reaches a stationary 
state. Quantities of interest are then averaged over an- 
other 1000 MC-steps, where one MC-step refers to the si- 
multaneous (parallel) update event of the state variables 
of the nodes. The used transition probabilities p^^' from 
state a to state a' in a single time step are Peh = 0.012, 
Phe = 0.05, Ppe = 0.25 and php is varied in the range from 
0.02 to 0.2 to produce the variation in /i = Php/Ppe- This 
procedure was repeated N times for different realizations 
of the graphs with N varying from = 50 for L = 2^"^ 
to A = 5 for L = 2^1. 

Fig. 121 shows how pp decays as a function of a host's 
parazitation probability parameter p. Below a size de- 
pendent critical value Pc{L) the parasites become ex- 
tinct resulting in a left-alone host population obeying 
dynamics defined by the SIS-model. For instance when 
L = 2^^ one may estimate that pc ~ 0.26. The inset 
in the Figure O strongly suggests that the relationship 
flj, ~ exp(— const / p) is established as in the SIS model 

To track Pc{L) more accurately we have studied the 
extinction probability VextiPi L) of the parasites during 



0.025 




0.08 0.16 0.24 0.32 0.4 0.48 0.56 0.64 



FIG. 3: Average parasite prevalence as a function of its 
spreading parameter p. The inset corroborates an Arrhenius 
relation pp ~ exp{— const / p) , as in the SIS model [13 . The 
error bars are smaller than the symbol size. 




0.07 0.08 0.09 0.1 0.11 

1/log(L) 



FIG. 4: Scaling of the critical point vs. system size. The 
dashed line works as a guide to the eye and suggest p.c{L) ~ 
l/log(L) as for the SIS model 

2000 MC-steps from different realizations of BA-graphs. 
The critical point is then determined to be the highest 
value of p below which the population dies away in a typ- 
ical realization of a BA-graph, and the sizes of the error 
bars in pc are estimated from the width of the window in 
which Vextip, L) decays from 1 to 0. Fig.^shows a scal- 
ing pc{L) ~ l/log(L) in the region 2^^ > L > 2^^, which 
again compares to the finite size scaling of the critical 
threshold in the SIS model |30l |. 

Since the probability for a node to become infected 
depends on its degree we next take a look at the para- 
site prevalence of nodes of degree k p^ in Fig. and the 
average degree of a site occupied by a parasite {k\p) in 
Fig.ini 

Fig. [51 shows that when approaching pc the relation- 
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Deqree (k) 



more and more difficult. Fig. |2| shows a consequence of 
this: the parasites do not prefer Uving in nodes of smaU 
degree anymore but, instead, the average degree of the 
nodes inhabited by them increases. In fact, as the inset 
of Fig.Elshows, the scahng {k\p) ^ 1/^ is estabhshed. A 
similar result should actually also hold for the SIS-model, 
and is predicted even by the mean-field equations H32|l 
and (|33|l . This in turn follows from the fact that for a de- 
creasing a the parasite density begins to saturate only at 
a higher and higher k (recall that it is linear in k for small 
degrees). We have also considered the autocorrelations 
of the time-series of the parasite prevalence, in analogy 
to Ref. [s^ . This decays exponentially with a time-scale 
constant that increases as the (pseudo-)critical point (for 
L fixed) is approached from above. 



FIG. 5: Average parasite prevalence and its dependence on 
the nodes degree. Here L = 2^® and only the pk of degree up 
to k — 100 are shown since the statistics for larger k become 
worse. 




0.48 0.56 0.64 



FIG. 6: The expectation value of the degree of a site occupied 
by the parasites. A scaling form for the average degree of 
parasitized nodes, {k\p), is found for small fi. The straight 
line in the inset is a guide to the eye. 



ship Pp ^ k begins to hold better and better whereas 
does not change remarkably since the host population is 
large. In fact, we have noted that, in analogy with the 
SIS-model, the scaling of Pp is not just a matter of coinci- 
dence but reflects the more general presence of the factor 
Pp ~ 1/(1 + const/fc) which is proportional to k at small 
p, or for large values of the constant. Generally, this be- 
havior implies that the largest connected component of 
hosts serves as a "scalefree" graph for the parasites thus 
partly explaining the absence of a critical point in the 
thermodynamic limit. 

As p —>■ pc, survival of the parasite population becomes 



V. DISCUSSION 

In this paper we have studied a two-population model 
("hosts" and "parasites"). First, as a preliminary, this 
problem was considered on the Bethe lattice. It turns out 
that the mean-field treatment can be augmented with the 
pair approximation. In particular we have been able to 
establish the generic form of phase diagram depicted in 
Fig- El This includes no tricritical point. 

The main finding of this work concerns the epidemic 
threshold of the parasites, in the presence of a non- 
zero host population, on scale- free graphs. Analytical 
arguments based on the neighbor-pair probabilities re- 
veal that in full analogy with the SIS model itself, the 
threshold is zero in the thermodynamic limit. Numeri- 
cal simulations on Barabasi- Albert model graphs imply 
that the finite-size behavior follows, also, the same scal- 
ing, and confirm this picture. These both findings might 
be surprising at first sight, due to the possible compli- 
cations from correlations. Concomitantly, correlations in 
the parasite dynamics are expected to follow the same 
picture as in the case of the SIS model. 

A striking feature related to correlated activity is the 
"escape" of the parasite population to vertices with, on 
the average, a high degree, which can actually be ex- 
plained within the standard picture of (SIS-type) popu- 
lation behavior as the prevalence is reduced by changing 
a control parameter. Due to the non-regular nature of 
the scale-free graphs we have not seen any indications of 
e.g. periodic, or chaotic oscillations that arise in many 
similar models on regular lattices p^ . [TTI | . Another pos- 
sible angle would be to study contact process-like models 
1^], where the spreading rate out of a graph vertex to a 
neighbor would depend on the degree of the out- vertex, 
for both parasites and hosts. The phase diagram of such 
model would be the same for the Bethe case, but for 
scale-free network one would, in analogy for the contact 
process itself |3l| , expect a finite threshold instead of the 
vanishing one for the SIS model. We have confirmed this, 
analytically, but obviously numerical studies would be of 
interest. 
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The results have imphcations, less for Bethe lattices 
which serve as an analytically tractable special case, but 
possibly for dynamical processes on real scale- free graphs. 
Examples can be found from ecology (metapopulation 
dynamics), where similar multi-species scenarios have al- 
ready been studied. Parasitoids do play a crucial role 
for the population dynamics of the endangered butterfly 
species melitaea cinxia in its fragmented habitat on the 
Aland islands in the Baltic Sea, which fit less well to sin- 
gle species models 1^- Due to the distribution of 
patch sizes and distances between them, the correspond- 
ing network model has a large tail degree distribution 
[33 ■ Whether a given patch is populated by hosts only 
or also by parasitoids depends on its local connected- 
ness. At least qualitatively the observations agree with 
those in Fig. Eland more systematic studies can be envi- 
sioned. The spreading rate depends on distance between 
and sizes of patches, in a nontrivial way |32| . If one trans- 
lates the underlying landscape to a network model, the 
resulting spreading rates may depend strongly or weakly 
on the degree of the emitting node, i.e. lie somewhere 
in the range between a generalized contact process and 
SIS-type models. Thus the limit considered by |3 1| may 
well be relevant in certain ecological systems. 

Another field of examples is epidemiology and vaccina- 
tion strategies. Knowledge on non-trivial network struc- 
tures in disease transmission can be used for vaccina- 
tion (see e.g. js^) or outbreak prediction, e.g. [sj, and 
also the importance of superinfections has been docu- 
mented (see e.g. a seminal work in an evolutionary con- 
text 38]). Our ansatz is an attempt to combine both 
points of view. From the scale-free network viewpoint 
the fundamental idea of concentrating the effort on nodes 
with a high k is valid here as well [3^ consider 



in particular the "escape" of parasites close to extinc- 
tion mentioned above. To fight parasites one needs, as 
well, to avoid random immunization. In this context an- 
other paradigmatic model is the "Susceptible-Infected- 
Removed" (SIR) model which is a variant of ordinary 
percolation. By taking in the HP model the right com- 
bination of limits for the parameters (essentially, disal- 
lowing recovery to the empty state from the H and P 
states) , one obtains a variant of the SIR which resembles 
in such language "bootstrap percolation" since the "R" 
( "P" ) sites are created only via contact with a neighbor 
in "R" . One should thus take note of possible general- 
izations of the HP model using similar recipes as can be 
applied to the SIR-style ones Hy. 

In the case of the SIS model, the cross-over (or the 
time-dependent picture) to the steady-state turns out to 
be interesting, which might be worth looking at here as 
well 01 ■ Another practical case related to this might 
be, say, viruses spreading as attachments to emails on 
the Internet [43 |. where again one is confronted with 
a dynamical graph (of email connections) on top of a 
larger one (internet). Finally, we would like to point 
out that our work could be extended to other similar 
multi-species models. An example would be a hierarchy 
of contact processes (A->B, B^C.) mill. 
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